/*
 *  VortSystem.h
 *  VortPan
 *
 *  Created by Gregor Veble on 7/14/10.
 *  Copyright 2010  . All rights reserved.
 *
 */

#ifndef VORTSYSTEM_H
#define VORTSYSTEM_H


#include <fstream>
#include <string>
#include <set>
#include "VortMesh.h"

#ifdef __APPLE__

#include <Accelerate/Accelerate.h>
#define pCPLX __CLPK_doublecomplex*
#define lapackint __CLPK_integer

#else

#define pCPLX Complex*


extern "C" {
	
	int dgetrf_(long* M, long* N, double *A, long *LDA, long* IPIV, long *INFO);
	
	int dgetrs_(char *TRANS, long* N, long* NRHS, double *A, long *LDA, long *IPIV, double *B, long *LDB, long *INFO );
}

#endif


namespace VortPan {
	
	struct ReflectedPos {
		Vector3D Position;
		Matrix3D Orientation;
		bool Parity;
	};
	
	
	typedef std::vector<re> Memory;
	typedef std::vector<Vector3D> MemoryVector;
	typedef std::vector<lapackint> IntMemory;

	typedef std::vector<ReflectedPos> PositionVector;
	
	class VortSystem : public VortBase {
	public:
		VortSystem(VortMesh& InitMesh);
		~VortSystem();
		
		VortMesh& GetMesh();
		
		void AddSymmetry(VortSymmetryPlane Plane);
				
		void SetVelocity(Vector3D Vel);
		
		re Potential(Vector3D Pos);
		Vector3D Velocity(Vector3D Pos, re Eps=1e-5);
		
		int CheckMesh(re Eps=0);
		
		int SolveSystem();
		Memory& GetDipolesReference();
		
		void StreamLine(Vector3D Initial, re timestep, re time, std::vector<Vector3D>& Points);
		void StreamLineRegular(Vector3D Initial, re timestep, re time, std::vector<Vector3D>& Points, re CoreSize);
		
		int ExportData(std::string FName, Memory& Data);
		
		int ExportClC(std::string FName, Vector3D Up=e3());


		int ExportIsWake(std::string FName);
		
		
		int ExportVectors(std::string FName, std::vector<Vector3D>& Data);
				
		void CalculateDirectPressures(Memory& Pressures);
		
		void CalculateDerivativePressures(Memory& Pressures, MemoryVector& Velocities, Vector3D* pForce=NULL, Vector3D* pMoment=NULL);
		
		void CalculateDerivativeVelocities(MemoryVector& Velocities);
		
		void CalculateAbsolutePressureGradient(Memory& Pressures, Memory& Gradients);

		void CalculateDirectedAbsolutePressureGradient(Memory& Pressures, MemoryVector& Velocities, Memory& Gradients);
		
		int RelaxWake(re Factor=ONE, re CoreSize=1e-3);
		
	protected:
		VortMesh Mesh;
		bool Symmetry;
		SymmetryPlaneVector SymmetryPlanes;

		void ReflectPositions(Vector3D Original, PositionVector& Reflections, PositionVector& Temp); //Reflections also contains the original
		
		re GetSymmetrizedSolidAngleInfluence(VortPolygon& Polygon, const Vector3D& Position, bool SelfInteraction);
		
		Vector3D GetSymmetrizedVelocityLatticeInfluence(VortPolygon& Polygon, const Vector3D& Position, re InvCoreSize2);
		
		re PolyToPolyMatrixElement(VortPolygon& PolySource, VortPolygon& PolyTarget);
		
		re WakeToPolyMatrixElement(VortPolygon& PolySource, VortPolygon& PolyTarget);

		re PolyToWakeMatrixElement(VortPolygon& PolySource, VortPolygon& PolyTarget, bool Sign);

		re WakeToWakeMatrixElement(VortPolygon& PolySource, VortPolygon& PolyTarget);
				
		void BuildMatrix();
		
		void BuildRHS();
		
		void DecomposeMatrix();
		
		void Solve(Memory& Vec);
		
		re FreePotential(Vector3D Pos);
		Vector3D FreeVelocity(Vector3D Pos);		
		
		re Potential(Vector3D Pos, Memory& Dipoles);
		Vector3D Velocity(Vector3D Pos, Memory& Dipoles, re Eps=1e-5);
		Vector3D VelocityAtSingularPoint(Vector3D Pos, Memory& Dipoles, re Eps=1e-5);
		Vector3D VelocityFromVortexLattice(Vector3D Pos, Memory& Dipoles, re CoreSize=1e-5);
		
		int RelaxWakePolys(VortMesh& OtherMesh, Vector3D Direction, re Factor, re CoreSize);
	
		int AssignSpace();
		int MInd(int i, int j);
		
		void ExportMatrix(std::string Name);
		
		
		bool MeshChecked;
		
		Memory Matrix;
		int SystemSize;
		Memory RHS;
		Memory DipoleStrength;
		IntMemory Pivot;
		
		Vector3D FarVelocity;
	};
	
}


#endif
